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PACS 05 . 20 . -y - Classical statistical mechanics 

PACS 61.20.Lc - Structure of liquids: Time-dependent properties; relaxation 

Abstract. - Slender rods in concentrated suspensions constitute strongly interacting systems 
with rich dynamics: transport slows down drastically and the anisotropy of the motion becomes 
arbitrarily large. We develop a mesoscopic description of the dynamics down to the length scale of 
the interparticle distance. Our theory is based on the exact solution of the Smoluchowski-Perrin 
equation; it is in quantitative agreement with extensive Brownian dynamics simulations in the 
dense regime. In particular, we show that the tube confinement is characterised by a power law 
decay of the intermediate scattering function with exponent 1/2. 



Brownian motion of highly anisotropic particles is con- 
siderably more complex than the diffusion of spherical 
objects, the basic understanding of which is founded on 
the seminal works by Einstein and von Smoluchowski. 
A shape anisotropy results in diffusion coefficients that 
depend on the direction of motion in the body frame, 
thus inducing a coupling of translation to the orienta- 
tion. This anisotropic dynamics has been investigated in 
recent experiments measuring diffusion coefficients of mi- 
^ crometre sized ellipsoids and rods by single particle track- 
^>ifig [1-3]; in particular, non-Gaussian statistics has been 
^observed [1]. Likewise in dynamic light scattering, the 
5^ Rotational and translational diffusion coefficients were de- 
termined simultaneously [4]. For these dilute systems, the 
ratio of diffusion parallel and perpendicular to the long 
symmetry axis was limited to values up to D»/D± « 4 in 
quasi two-dimensional (2D) confinement. 

Considerably higher values of this ratio have been ob- 
served in simulations of semi-dilute suspensions of slender 
rods, yielding D\\/D±_ up to values of 50 [5,6]. This in- 
crease in anisotropy is caused by the steric constraints 
imposed by surrounding rods; thereby the transverse and 
rotational motion is suppressed, whereas the longitudinal 
transport is barely influenced [7] . An intermediate regime 
of anisotropic diffusion has been derived for ballistic nee- 
dles within kinetic theory [8] and was observed in simula- 
tions [9]. 

For a finite width b, the rods undergo a phase tran- 
sition to the nematic phase at densities of the order 



of 1/bL 2 as predicted by Onsager [10]. The dynamics 
in this ordered phase becomes trivially anisotropic and 
splits into a fast diffusion along the nematic director axis 
and slower diffusion perpendicular to it. Such a pro- 
nounced anisotropic diffusion has been observed in simula- 
tions of nematic elongated ellipsoids [11] and spherocylin- 
ders [12,13]. Experiments have also clearly demonstrated 
orientation-dependent diffusion in colloidal nanorods in 
the isotropic and nematic phase [14] and in various liquid 
crystalline phases of fd viruses [15,16]. The phenomena 
connected with the nematic phase transition are essen- 
tially understood and will not be discussed in this work. 

In the isotropic phase, experiments and computer simu- 
lations have been restricted to determine the diffusion co- 
efficients from the mean-square displacements. The com- 
plex interplay of translational and rotational motion as ex- 
emplified for a single free ellipsoid [1] has not been studied 
for the strongly hindered motion in solution yet. In prin- 
ciple, one should characterise the dynamics in terms of 
an intermediate scattering function or a van Hove corre- 
lation function, as has been done recently for the smectic 
phase [16]. 

A general theory for the anisotropic motion of rods in 
entangled suspensions is a long-standing problem, due to 
the intricacy of the many-body interaction. In particu- 
lar, such a theory should include memory effects from the 
translation-rotation coupling and explain the emergence of 
new macroscopic time and length scales. To account for 
phenomena which depend on the simultaneous interaction 
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with many particles, a non-perturbative approach is re- 
quired. Substantial progress would be achieved in terms 
of an effective one-particle theory that allows for quanti- 
tative predictions. 

Unhindered motion. — Already at the level of a sin- 
gle, free rod the problem is involved: the dynamics of 
the probability distribution ^(i?, m, t) is governed by the 
Smoluchowski-Perrin (SP) equation [17-19], 

dt^ = -D rot LH + 8 R • [(£>|| - D ± )uu + D±l] ■ d R ^, 

(1) 

with the centre of mass position of the rod R, its orien- 
tational unit vector u, the rotational diffusion coefficient 
D rotl and the angular part of the Laplacian, — L 2 . The 
full formal solution of this equation in 3D was given in 
ref. [20]; its quantitative evaluation is still missing, al- 
though certain aspects are well understood [18,19]. To 
the best of our knowledge, the solution ^>(R,u,t) of the 
Smoluchowski-Perrin equation has not been discussed for 
2D systems. Equation (1) is trivially solved for isotropic 
diffusion, D± = Du, and at macroscopic time scales, 
t> 1/ D rot , where the translation-rotation coupling is re- 
laxed. Then, the second term may be averaged over u, 
yielding the average diffusion coefficient. 

Model for interacting rods. — We study a model 
for the sterically constrained dynamics of rods that natu- 
rally induces anisotropic diffusion of arbitrarily large ra- 
tios D\\/D±. The model is set up in a 2D environment; 
the dynamic properties we focus on are however equally 
present in a 3D embedding space. We demonstrate that 
the Smoluchowski-Perrin equation provides an excellent 
effective theory for the dynamics, when the measured dif- 
fusion coefficients serve as input parameters. In particu- 
lar, we compare mean-square displacements (MSDs) and 
the intermediate scattering function from computer sim- 
ulations with their exact results from eq. (1). Further- 
more, we have discovered an intermediate algebraic decay 
in the intermediate scattering function, characteristic for 
the anisotropic sliding motion. 

The model considers the overdamped motion of a sin- 
gle rod with zero width exploring a plane with randomly 
distributed, hard point obstacles [21]. Then the orienta- 
tional unit vector of the rod is parametrised by one angle, 
u(t) = (cos$(i),sin$(i)), and L 2 reduces to — 9J,. The 
model ignores any excluded volume to exclusively concen- 
trate on entanglement effects. The physical properties of 
this system are thus controlled by a single parameter, the 
reduced density n* := riL? , where L denotes the length 
of the rod, and n the number density of the obstacles. 
The model is closely related to 3D suspensions of rods 
when considering a planar section. A tagged rod found 
initially in this plane is approximately confined to it by 
the neighbouring rods for times shorter than the orienta- 
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tional relaxation time, t <C T rot 
scale present in the system [21]. 



= D TOt , the largest time 
The constrained motion 



Fig. 1: (colour online) Simulated body- frame MSD in the en- 
tangled regime. For each density, three observables are shown: 
The MSD parallel to the rod's axis (topmost, dash-dotted), 
perpendicular (solid), and the MS AD (broken). Arrows at the 
bottom indicate T con f for the corresponding density; the diffu- 
sive regime is shaded. 



of the rod then corresponds to diffusion in a planar course 
of localised intersection points. 

Anisotropic diffusion. — The microscopic motion 
of the rod is diffusive, with diffusion coefficients chosen 
according to first order hydrodynamics of a slender rod, 
D° ± = D{|/2, and D r ° ot = 6L>j|/£ 2 . For the computer sim- 
ulations we have combined the Langevin equations cor- 
responding to cq. (1) with an event-driven algorithm to 
detect the collisions between rod and obstacles [21-23]. 
The MSDs in fig. 1 visualise essential properties of the 
model in the semi-dilute regime, n* 3> 1. We define 
the displacement in the body-fixed frame along the axis as 
Ai?||(i) := J Q R(t') ■ u(t')dt' , and similarly the transverse 
part AR±(t). The parallel MSD, Srf(t) := (AR^t) 2 ), is 
not affected by the obstacles at all, due to zero excluded 
volume. Consequently, the parallel diffusion coefficient is 
independent of the density, D\\ = DjJ. In contrast, the 
perpendicular MSD 6rj_(t) and the mean-square angular 
displacement (MSAD) <5i9 2 (i ) enter a plateau beyond a 
density-dependent time scale T con f. The plateau reflects 
the local confinement to an effective cage built up by the 
surrounding obstacles, referred to as "tube" [19]. Its di- 
ameter d is determined by nLd 1, leading to a rela- 
tion for the time when the confinement becomes effective, 
r conf := d 2 /D\ w l/n 2 L 2 D° ± . At the time scale r = 
L 2 /D||, the rod moves a distance comparable to its length 
L, hence it leaves the tube and the MSDs become diffusive 
again. From the long-time asymptotes, the diffusion co- 
efficients are read off, e.g., D±(n*) = lim t _ 00 d t 5r 2 _(t)/2. 
Fig. 2 demonstrates the huge suppression of perpendicu- 
lar and rotational diffusion coefficients: both scale with 
obstacle density as n~ 2 , as has been argued by Sza- 
mel [7]. A consequence of the vanishing perpendicular 
component is the saturation of the centre of mass diffusion, 
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Fig. 2: (colour online) Density dependence of the diffusion coef- 
ficients. Left axis: rotationai diffusion and transverse diffusion; 
right axis: centre of mass diffusion. 



D cm {n* > 1) = 2L>" m /3. For the largest simulated den- 
sity, n* = 2240, the anisotropy ratio approaches a value 
of D\\/D±_ ~ 10 5 . Note that these large ratios in dense 
systems do not require an anisotropy on the microscale; 
rather the anisotropic motion is generated dynamically 
from the strong interaction with many obstacles. 

A characteristic crossover between anisotropic and 
isotropic dynamics is seen clearly when plotting the time- 
dependent diffusion coefficients in a space-fixed frame, 
D x (t) := 9 t (Ax(t) 2 ) f /2 and D y {t) := d t (Ay{t) 2 ) { /2, 
where the subscript 'f indicates that the initial orienta- 
tion is fixed to the x-axis; seen fig. 3. The memory of 
the orientation is lost only at times larger than r ro t ~ ?i 2 , 
resulting in a time window of anisotropic diffusion that 
is significantly prolonged with increasing density. Note 
that fig. 2 directly visualises the density dependence of 
the timescale indicating the crossover to isotropic dynam- 



ics, T rot = Arof 

For comparison with the unhindered anisotropic mo- 
tion, eq. (1) is solved for the conditional probabil- 
ity distribution ^(R, £|$o); with the initial condition 
ty(R,ti,t = 0|i? ) = 6(R)8(d - j?o)- A Fourier trans- 
form defines the characteristic function Gfc(i?,i|i?o) = 
fe~ lk R $>(R,'&,t\'&o)dR. Its equation of motion attains 
the form of a Schrodinger equation, 



dtGk — —HoGk — VGk, 



(2) 



with the operators 7io = — D Ioi d^ and V = (D\\ — 
D±)(k ■ u) 2 + D±k 2 . Perturbation theory in kL <C 1 now 
solves iteratively for Gk- The central quantity of inter- 
est is the intermediate scattering function, F(k, i|$o) : = 
/G fe (^,i|i?o)di?. Up to fourth order in k, 



F(fc,t|#o) = l-DkH-^n(t) 



Fig. 3: (colour online) Deviations of the diffusion coefficients 
D x (t) and D y (t) in the space-fixed frame from the isotropic 
value D. The initial orientation is fixed to the x-axis; normal- 
isation is chosen such that values of ±1 indicate diffusion with 
Du and D±, respectively. Symbols show results of simulations, 
solid lines the effective Perrin theory; arrows indicate r ro t for 
the different densities. The broken line displays the measured 
[D±(t) — D]/D a at n* — 12 for comparison. For t — > 0, the 
perpendicular diffusion in the simulations asymptotically ap- 
proaches [D y (t -> 0) - D]/D a = [1 - D l[ /D ± (n*)}- 1 . 
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The notation is abbreviated with the isotropic and 
anisotropic diffusion coefficients, D = (D\\ + D±)/2, 
D a = (_D|| — D±)/2, the wavevector components k± = 
(k x ±ik y )e^°, and Tj(t) := J *e^ D ™ tS ds. Note that 
the dependence of eq. (3) on $o is hidden in the defini- 
tion of k± — this property originates from rotational sym- 
metry and holds to all orders in k. From F(k,t\-do) all 
moments are obtained by derivatives, e.g., (Ax(t) 2 ){ = 
— dl^F(k, i|#o)U=o- Some of these moments have been 
calculated by Han ct al. [1] from the Langevin equations 
equivalent to eq. (1). 

Effective Perrin theory. — The idea is to use the so- 
lutions of the Smoluchowski-Perrin equation with the dif- 
fusion coefficients Dm, D±(n*), D Iot (n*) — measured in the 
simulations — to obtain a prediction for time-dependent 
MSDs and the intermediate scattering function in the pres- 
ence of obstacles. In other words, the results from the 
simulations for two transport coefficients in semi-dilute 
systems are employed in the theory derived for free dif- 
fusion to describe the full time-dependence of a semi- 
dilute suspension. Since the orientation changes only 
gradually, the long-time rotational motion is described 
by diffusion [21]. Thus the description in terms of an 
effective Perrin theory has to be trivially correct in the 
relaxed regime where translational diffusion is isotropic, 
i.e., on macroscopic time and length scales, t > r ro t and 
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k 1 > irot := V-DaTrot- 

We will show that in fact the effective Perrin the- 
ory constitutes a quantitative mesoscopic theory in the 
dense regime, n* 3> 1, and successfully describes also 
the translation-rotation coupling induced by the interac- 
tion with the many obstacles. That such an approach 
should work in principle has been anticipated earlier for 
dense needle liquids [8]. A comparison with simulated 
MSDs. fig. 3, reveals excellent agreement down to the 
time scale tq. For short times, f C To, the space- and 
body-fixed frames coincide, implying D x (t) ~ D\\(t) and 

Dy{t) ~ D X (t). 

To access the full range of wavenumbers, we construct 
the exact solution of G k ($, f |$o) m terms of Mathieu func- 
tions. The equation of motion, eq. (2), is rewritten as 

d t G k = D mt dlG k ~k 2 (D + D a cos 2d) G k , (4) 

in a coordinate frame with fe = ke x . A separation 
ansatz, Gk($,t) = (?fc(i?)e~ At , yields the Mathieu equa- 
tion, = d 2 g k + (a — 2g cos 2^) g k , with the parameter 
q = D a k 2 /2D Iot and the eigenvalue a = (A — Dk 2 )/ D TOt . 
The general solution is thus a linear combination of even 
and odd eigenfunctions, cej(-d,q) and sej(i?, q) [24]; 1 the 
decay rate A = A(a, k) = aD rot + Dk 2 depends on the 
corresponding even and odd eigenvalues, a — > aj(q) and 
a — > bj(q), respectively: 



G k (0,t\0 o ) 



J2 [e- A ^' fe ) t ce J -(t? ,g)ce J -(t?, g ) 



- x ^^ t se j ^ ,q)se j (^q)\ . (5) 



The intermediate scattering function with unconstrained 
initial orientation is obtained by integrating over ■& and 
0o, 



F(k,t) 



,-k 2 Dt c -a 2j (q)D Iat t 



A^\q) 



(6) 



3=0 



with coefficients A^\q) := J ' 27r cc2j(i9, q)&$/\/2~K. 

The convergence of the sum is determined by the 
magnitude of q. The eigenvalues a,j{q) are ordered 
ascendingly in j; furthermore, the coefficients fulfil 

^0 (?) = C (?"')■ Hence for q <C 1, the low-j terms 
yield the major contributions to the sum, F(k,t) = 



-k'Dt 



rep- 



c q ZW/2(J _ ? 2/ 8 ) + c -4D rottg2/ / 8 + 0^4) 

resenting the first correction to isotropic diffusion. Thus 
q 2 < 1 defines the macroscopic regime, corresponding to 
(/cL rot ) 4 < 4, and L ro t is derived as the relevant macro- 
scopic scale. 

The opposite limit, q 1, is relevant for large densities 
due to the suppression of D TOt . The asymptotic expan- 
sion a,j(q) ~ — 2q + (4j + 2)^fq + 0(1) reveals a reduc- 
tion of the exponential prefactor in cq. (6) to c~ k D±t , 

1 We use the normalisation j 2lT ce^ q)d-& = 1 and 



and a large number of terms contribute. Then, the termi- 
nal relaxation is ruled by an exponential with decay rate 
r t -J m := kV2D a D mt + k 2 D ± . 

Based on the simulation results for the intermediate 
scattering function, we test the range of validity of the 
effective Perrin theory in fig. 4. The quality of our data, 
approaching a signal-to-noise ratio of 10 -4 in the scat- 
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Fig. 4: (colour online) Time and wavenumber dependence of 
the intermediate scattering function for densities (a) n* = 0, 
(b) n* = 4, (c) n* = 40, and (d), (e) n* = 400. Symbols 
represent simulation results, lines the effective Perrin theory, 
eq. (6), and time is in units of k 2 D. By rescaling time with 
Tterm, inset (e) visualises the terminal relaxation. 
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Fig. 5: (colour online) intermediate scattering function for fixed 
wavenumber and varying density. Increasing confinement is 
manifested in the simulations (symbols) by the development 
of a power law decay. Lines show the effective Perrin theory; 
the perturbative correction to eq. (7) is indicated by the thick 
black line for n* = 1200. 

tcring function, allows for a clear distinction of features 
on a large range of time and length scales. 2 By con- 
struction, the theory describes the data at zero density 
[panel (a)]. For intermediate densities, n* ~ 1, the Per- 
rin prediction strongly deviates from the simulation. The 
two lowest wave numbers are in the macroscopic regime; 
correspondingly the curves coincide after proper rescal- 
ing [panel (b)], and the Perrin prediction applies in the 
trivial sense. At larger densities, the agreement becomes 
increasingly accurate as the wave number decreases [pan- 
els (c), and (d)] also for wavenumbers beyond the macro- 
scopic regime. Second, deviations are shifted to larger 
wavenumbers for higher densities. We conclude that there 
is a density-dependent length scale £ determining the va- 
lidity of our coarse-grained approach. From all simulated 
densities, we have identified £ as the typical distance be- 
tween obstacles, £ := -nT 1 ! 2 : for the shown densities, 
n* = {4,40,400}, the corresponding wavenumbers 27r/£ 
are 13L -1 , 40L -1 and 126L -1 , respectively The slow- 
ing down of the dynamics with increasing confinement is 
displayed in fig. 5. For high wavenumbers k ~ 27r/d the 
intermediate scattering function probes the formation of 
the tube at time scales T con f, sec figs. 4c and 5. Once the 
tube confinement becomes effective, an intermediate alge- 
braic decay emerges, F(k,t) ~ i -1 / 2 , which we attribute 
to the sliding motion inside the tube. This power law is 
cut off by an exponential relaxation at Tterm, see fig. 4e. 

Origin of the power law. — The power law is hidden 
in eq. (6) in the sum of many exponentials for <; > 1. For 
strongly suppressed perpendicular and rotational motion, 

2 To achieve sufficient statistics, we collected at least 350 trajec- 
tories for each density. For the largest densities, the simulation of 
a single trajectory took about 13 days of CPU time on a AMD 
Opteron® 2.6 GHz core. 



eq. (4) is approximated by d t Gk = —k 2 D^cos 2 (d) Gk, 
which yields 

F(k,t)=e- k2D ^ 2 I (k 2 D ll t/2) ~ (7rfc 2 ^||t)" 1/2 . (7) 

The second relation results from an expansion of the mod- 
ified Bessel function of the first kind Iq(z) for large argu- 
ment; it quantitatively reproduces the scattering function 
in the power law regime, shown in fig. 4d. Perturbation 
theory in D rot also captures the terminal relaxation with 
the previously calculated Tterm, see fig. 5. 

Conclusions. — From our analysis we conclude that 
the motion of thin rods in concentrated suspensions ex- 
hibits a rich interplay of time and length scales, as ex- 
emplified in the MSDs and intermediate scattering func- 
tion. The hindered motion leads to a strong anisotropic 
dynamics manifested in a significant translation-rotation 
coupling which persists up to macroscopic length and time 
scales, L rot and r rot . This dynamically induced coupling 
is an emergent phenomenon with long memory effects. 
The confinement in the tube gives rise to a well sepa- 
rated spectrum of length scales d -C £ -C L -C £ ro t- 
The tube becomes effective once the rod encounters new 
steric constraints — thus the interparticle distance £ consti- 
tutes the lower length scale of a mesoscopic window, where 
an effective theory with renormalized parameters becomes 
valid; this regime extends up to L mt . Note that such a 
window opens only in the strongly anisotropic regime and 
is absent for, e.g., ellipsoids with moderate aspect ratio. 
There the orientational motion is dominated by excluded 
volume effects, and becomes neither diffusive nor expo- 
nential [25,26]. For strongly entangled suspensions, the 
orientation changes only gradually by tube renewals, and 
it has to be included in the set of slow degrees of freedom 
in addition to the translation. The effective Perrin the- 
ory constitutes a Markov process in these variables, and 
the long memory observed in the intermediate scattering 
function is generated by integrating out the slowly vary- 
ing orientation. This mesoscopic description is indepen- 
dent of the details of the tube generation; in particular, 
the obstacles may also fluctuate in time and space, or even 
disappear. 

The Perrin approach fails to capture the fact that 
strongly confined needles have to diffuse along their axis a 
distance L to relax the tube constraint. However, since the 
rotational diffusion D mt is much slower than D\\/L 2 , this 
appears to be negligible, at least on the scales investigated 
here. 

We expect that our findings are also relevant for slen- 
der rods of finite width w. Then the motion of a sin- 
gle rod is again strongly confined by a tube comprised 
of the surrounding particles. For large aspect ratios L/w 
the isotropic-nematic spinodal occurs at 3D concentrations 
c w A/wL 2 [22], e.g., at c = 200/L 3 for L/w = 50. For 
densities 100 < cL 3 < 150, 3D simulations of suspen- 
sions of rods with this aspect ratio exhibit significantly 
anisotropic diffusion of about 25 < Dn/D± < 50 [5,6]. 
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In our 2D model such ratios are attained for n* « 20, 
already at the onset of the regime of entangled diffusion. 
Experimental realisations of highly entangled rod suspen- 
sions may be achieved, e.g., for the tobacco mosaic virus. 
Mutants of this rod-shaped virus have been observed with 
an aspect ratio of 50 and higher [27] , thus entangled sus- 
pensions below the nematic phase transition should be 
feasible. Another promising model system are solutions 
of microtubules; their aspect ratio can even be consider- 
ably larger [28, 29], and with a ratio of 100 they can still 
be considered as approximately stiff rods. A third exam- 
ple are solutions of fd viruses, exhibiting a high degree 
of monodispcrsity and a comparable aspect ratio of about 
130 [16]. 

In the preceding discussion wc have neglected hydrody- 
namic interactions mediated by the solvent. Inclusion of 
these effects present a considerable challenge both to the 
theory and the simulation techniques, and is beyond the 
scope of our work. The possible influences of hydrody- 
namic forces are twofold: First, a logarithmic length de- 
pendence, ln(L/w), decorates the microscopic expression 
for the diffusion coefficients. Since the ratios of D±, Du, 
and D mt , arc independent of this correction, this effect 
is already accounted for by choosing To = L 2 / 'D?, as the 
basic unit of time, as we have done here. Second, nonaddi- 
tive hydrodynamic forces between the rods might change 
the observed dynamics, in particular at short times [30]. 
Recently, Pryamitsyn and Ganesan [31] have argued on 
the basis of computer simulations that effects of hydro- 
dynamic interactions on the translational and rotational 
diffusivity are secondary relative to the steric interactions. 
Assuming that hydrodynamic interactions merely renor- 
malise the macroscopic transport coefficients, the mecha- 
nism of translation-rotation coupling remains unaffected, 
and our theory should still be applicable. 

It is straightforward to extend the concept of an ef- 
fective Perrin theory to 3D suspensions of long, thin rods 
since the motion in the tube is essentially one-dimensional. 
The formal solution is then provided in terms of spheroidal 
wave functions [8,20], which serves as a starting point to 
calculate the intermediate scattering function. In partic- 
ular, the confined motion of rods is characterised again 
by a power law decay of the intermediate scattering func- 
tion, which we predict to F(k,t) ~ (4fc 2 D||i/7r) -1 / 2 . This 
algebraic decay constitutes a generic feature of the slid- 
ing motion and should be observable directly in scattering 
experiments. 

* * * 
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